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Particles in the atmosphere reflect incoming sunlight, tending to cool the Earth below. Some 
particles, such as soot, also absorb sunlight, which tens to warm the ambient atmosphere. 

Aerosol optical depth (AOD) is a measure of the amount of particulate matter in the atmosphere, 
and is a key input to computer models that simulate and predict Earth’s changing climate. The 
global AOD products from the Multi-angle Imaging SpectroRadiometer (MISR) and the 
MODerate resolution Imaging Spectroradiometer (MODIS), both of which fly on the NASA 
Earth Observing System’s Terra satellite, provide complementary views of the particles in the 
atmosphere. Whereas MODIS offers global coverage about four times as frequent as MISR, the 
multi-angle data makes it possible to separate the surface and atmospheric contributions to the 
observed top-of-atmosphere radiances, and also to more effectively discriminate particle type. 
Surface-based AERONET sun photometers retrieve AOD with smaller uncertainties than the 
satellite instruments, but only at a few fixed locations. So there are clear reasons to combine 
these data sets in a way that takes advantage of their respective strengths. 

This paper represents an effort at combining MISR, MODIS and AERONET AOD products over 
the continental US, using a common spatial statistical technique called kriging. The technique 
uses the correlation between the satellite data and the “ground-truth” sun photometer 
observations to assign uncertainty to the satellite data on a region-by-region basis. The larger 
fraction of the sun photometer variance that is duplicated by the satellite data, the higher the 
confidence assigned to the satellite data in that region. In the Western and Central US, MISR 
AOD correlation with AERONET are significantly higher than those with MODIS, likely due to 
bright surfaces in these regions, which pose greater challenges for the single-view MODIS 
retrievals. In the east, MODIS correlations are higher, due to more frequent sampling of the 
varying AOD. These results demonstrate how the MISR and MODIS aerosol products are 
complementary. The underlying technique also provides one method for combining these 
products in such a way that takes advantage of the strengths of each, in the places and times 
when they are maximal, and in addition, yields an estimate of the associated uncertainties in 
space and time. 
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10 ABSTRACT 

1 1 The Multi-angle Imaging SpectroRadiometer (MISR) and the Moderate Resolution Imaging 

12 6SHF WURUDGLRPHWHU □ □ 02',6 □ □ DERDUG □ WKH □ 1 $6$ n (DUWK n 2E VHUYD WLRQ n 6\VWHPtV n 7HUUI 

13 measuring aerosol optical thickness (AOT) since early 2000. These remote-sensing platforms complement 

14 the ground-based AErosol RObotic NETwork (AERONET) in better understanding the role of aerosols in 

15 climate and atmospheric chemistry. To date, however, there have been only limited attempts to exploit the 

16 complementary multiangle (MISR) and multispectral (MODIS) capabilities of these sensors along with 

17 the ground-based observations in an integrated analysis. This paper describes a geostatistical data fusion 

18 technique that can take advantage of the spatial autocorrelation of the AOT distribution, while making 

19 optimal use of all available datasets. Using Level 2.0 AERONET, MISR and MODIS AOT data for the 

20 contiguous US, we demonstrate that this approach can successfully incorporate information from multiple 

21 sensors, and provide accurate estimates of AOT with rigorous uncertainty bounds. Cross-validation 

22 results show that the resulting AOT product is closer to the ground-based AOT observations than either of 

23 the individual satellite measurements. 



24 


1 INTRODUCTION 


25 Atmospheric aerosols play an important and dynamic role in climate and atmospheric chemistry. The 

26 climatic effects of aerosols had already been recognized in the 1970s [Andreae, 1995] but the focus of 

27 scientific attention shifted only during the late 1980s due to the impact of the growing concentrations of 

28 C0 2 and other greenhouse gases. Although the radiative forcing of aerosols is still highly uncertain 

29 [IPCC, 2007], it is well understood that aerosols contribute significantly to reflected solar radiation (the 

30 aerosol direct effect) and modify cloud properties (the aerosol indirect effect), producing a net cooling of 

31 the Earth surface, and can also absorb sunlight, thereby warming the ambient atmosphere. Because 

32 aerosols have short atmospheric lifetimes of about a week [Andreae, 1986], they have a heterogeneous 

33 spatial and temporal distribution. Accurately capturing this heterogeneity, and assessing the impact of 

34 tropospheric aerosols on regional and global energy budgets, therefore requires diumally resolved 

35 observations from some combination of satellite and suborbital measurements. 

36 Two space-based instruments that aim to fulfill this need are the Multi-angle Imaging SpectroRadiometer 

37 (MISR) and Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the NASA Earth 

3 8 ObservDWLRQ □ 6WWHPKV □ 7HTJMJi0MiMdffi3(M J JMHve observations of the tropospheric aerosol 

39 optical thickness (AOT), among other parameters [Diner et al., 1998; Kaufman et al., 1997]. Column 

40 AOT is defined as the integral of aerosol extinction from the surface to the top of the atmosphere. 

41 Although these two sensors are on the same platform, discrepancies exist between them in retrieved AOT 

42 over both land and ocean regions [Penner et al., 2002; Myhre et al., 2005; Kinne et al., 2006]. These 

43 discrepancies are due to the differences in assumptions in the retrieval algorithms [Kahn et al., 2007], 

44 observed wavelengths and viewing geometries [IPCC, 2007], and the spatial resolution of observations 

45 [Xiao et al., 2009], among other reasons. Methods for evaluating data from these and other instruments 
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46 are needed, as are approaches for assessing the information content of these data for providing the best 

47 possible representation of the spatial and temporal variability in AOT. 

48 The common way of validating the satellite AOT retrievals has been through the use of the ground-based 

49 Aerosol Robotic Network (AERONET; Holben et al., 1998), which provides sparse but relatively reliable 

50 AOT observations. Comparisons between AOT retrieved from space-based instruments and AERONET 

5 1 data have been used in a variety of contexts to explore the similarities and differences between the MISR 

52 and MODIS products. These comparisons have focused on MISR and AERONET [Liu et al., 2004a; 

53 Kahn et al., 2005a, 2005b; Jiang et al., 2007; Chen et al.. 2008], or MODIS and AERONET [Chu et al., 

54 2002; Levy et al., 2003, 2005; Renter et al., 2005], and have been specifically targeted at refining the 

55 retrieval algorithms of the individual sensors for different aerosol regimes. 

56 Several studies have also looked at the discrepancies between MISR and MODIS [e.g., Abdou et al., 

57 2005; Liu et al., 2007; Prasad and Singh, 2007; Vermote et al., 2007; Xiao et al., 2009; Kahn et al., 

58 2009], mostly by comparing them with the AERONET measurements. These studies have concluded that 

59 the major differences can be attributed to location (for example, retrievals near aerosol source regions 

60 and/or presence of clouds, retrievals over land versus water) and the aerosol retrieval algorithms over 

61 those locations. Recently, Kahn et al. [2009] compared MISR and MODIS datasets, and found strong 

62 correlations of 0.9 and 0.7 between MISR and MODIS over ocean and land, respectively. Discrepancies 

63 between the instruments were traced back to sampling differences, known algorithmic issues, or other 

64 mechanisms contributing to aerosol retrieval error. Some of these mechanisms that have been highlighted 

65 previously are aerosol model differences [Abdou et al., 2005; Kahn et al., 2007], presence of clouds 

66 [Martonchik et al., 2004; Kahn et al., 2007; Xiao et al. 2009; Kahn et al., 2009], dust [Kalashnikova and 

67 Kahn, 2006; Martonchik et al., 2004], biomass burning [Kahn et al., 2005a; Chen et al., 2008], and other 
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68 biospheric and anthropogenic factors [Prasad and Singh, 2007; Xiao et al., 2009]. Statistical comparisons 

69 have also been carried out between MISR, MODIS and AERONET by Liu and Mishchenko [2008] and 

70 Mischenko et al. [2009], although some of the statistical techniques used have subsequently been 

71 questioned [e.g., Kahn et al., 2009]. Overall, the existing literature has resulted in a complex set of 

72 conclusions regarding the ways in which MISR, MODIS, and AERONET record AOT [Xiao et al., 2009]. 

73 For example, Liu et al. [2007] conclude that MODIS generally retrieves higher AOT relative to MISR 

74 over land, whereas both MODIS and MISR tend to underestimate AERONET AOT measurements for 

75 AOT higher than about 0.5. Similar underestimation is reported by Jiang et al. [2007] and Kahn et al. 

76 [2005a], whereas others conclude that MISR overestimates AERONET AOT observations over water 

77 [e.g., Abdou et al., 2005; Kahn et al., 2005b; Liu et al., 2004a]. 

78 Given the limitations inherent to each of the available data streams, combining multiple data types may 

79 provide an opportunity to optimally estimate the spatial and temporal distribution of AOT. Some studies 

80 have found the correlation between the AOT data from multiple sensors to be sufficiently strong to justify 

81 the use of ground-based and space-based observations together [Liu et ah, 2004a; Prasad and Singh, 

82 2007; Jiang et al., 2007]. However, most of the data fusion attempts have been limited to merging data 

83 from multiple space-based instruments, including Level IB (i.e. radiance) data [Loeb et al., 2006], Level 

84 2 data of geophysical parameters [Gupta et al., 2008] and aerosol optical depth [Nguyen 2009], and 

85 gridded level 3 datasets [Acker et al., 2007]. Recently, Kinne [2009] presented an approach for integrating 

86 a weighted composite of remote sensing AOT observations with AERONET AOT through an empirical 

87 averaging procedure. 

88 Given the complementary capabilities of the AERONET, MISR and MODIS sensors (see Section 2), it 

89 seems natural to investigate whether it is possible to merge data from these different sensors in a 


Page | 5 



90 statistically rigorous framework to obtain an improved AOT product. Such a product could be used to 

91 address scientific issues related to air quality and the radiative effects of aerosols, and in particular, be 

92 used to evaluate model predictions of aerosol distributions. 

93 The objective of this work is to investigate the applicability of Universal Kriging, a simple geostatistical 

94 data fusion approach, for merging multiple AOT datasets. The approach yields a statistical best-estimate 

95 of the AOT spatial distribution, together with a quantification of the associated uncertainty. The estimated 

96 AOT distribution is based only on the available AOT data, and does not incorporate information or 

97 assumptions about atmospheric transport or source regions. Given that the availability of multiple satellite 

98 datasets has already resulted in a research shift from modeling-only to observational-based assessments of 

99 aerosol forcing [Yu et al., 2006], geostatistical data fusion can potentially provide useful optimal fused 

100 datasets, taking advantage of the strengths, and minimizing the limitations, of each individual sensor in a 

101 new way. 

1 02 The remainder of this paper is organized as follows. Section 2 provides a description of the MISR, 

103 MODIS and AERONET data used in the presented analysis. Section 3 gives an overview of the applied 

104 method and examined test cases. Results are presented and discussed in Section 4. 

105 2 DATA 

106 The description of the datasets presented here covers only the specific data products used in this study. 

107 The reader is referred to Martonchik et al. [2009] and Remer et al. [2009] for descriptions of the retrieval 

108 algorithms and Yu et al. [2006] for an overview of how tropospheric aerosols are measured. All analyses 

1 09 are performed using data from 2001. 
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2.1 AERONET 


111 AERONET is a globally distributed network of over 200 automated ground-based instruments covering 

112 all major tropospheric aerosol regimes [Holben et al., 1998; 2001]. The instruments used are CIMEL 

113 sun/sky radiometers that make direct sun measurements with a 1 .28UFull field-of-view every 15 minutes 

1 14 in eight spectral bands [Holben et al., 1998]. Level 2 (validated) AOT data are used here for 32 sites 

115 within the continental United States. The AERONET data archive (http://aeronet.gsfc.nasa.gov), includes 

AOT at diffHUHQW □ ZDYHOHQJWKV □ □ UHODWLYH □ HUURUV □ RI □ $27 □ □ $Q JVWURP □ H[SRQHQ WV □ □ 



129 resolution of the dataset is 17.6 km. Theoretical sensitivity studies for MISR [Kahn et al., 2001] have 

130 estimated the standard deviations of the measurement error associated with the optical depth to be 0.05 

□RUD □ □ □ 



150 The AOT data from the three sensors cannot be compared directly, in part because they are reported at 

151 different spatial resolutions. Therefore, following the methodology of Liu et al. [2004a], the mean of 

1 52 MISR and MODIS observations within a 0.5° by 0.5° bounding box around each AERONET site is used 

153 as a basis for comparison to AERONET data, which are themselves averaged over ±30 min from the 

1 54 Terra overpass. Correlation coefficients are used to characterize the agreement between daily data-pairs 

155 from the 32 AERONET sites and the corresponding MISR or MODIS observations at those sites. 

156 3.2 Investigation of the Spatial and Temporal Variability in MISR and MODIS AOT 

1 57 AOT varies spatially and temporally. This variability can be quantified using variogram analysis, a 

158 geostatistical spatial analysis tool. Although the AERONET network is too sparse to independently 

159 characterize the spatial variability at the continental scale, it can be used for regional analyses in areas 

160 when the network is relatively dense. On the other hand, the dense MISR and MODIS data coverage 

161 provides good information about AOT spatial variability as captured by these instruments. Analysis of the 

162 MISR and MODIS AOT spatial variability provides insights into differences in the way that these 

163 instruments capture the AOT distribution. Differences may be due to the differences in the observational 

164 spatial resolution and sampling, instrument signal-to-noise ratios, or retrieval algorithms. 

165 For assessing the AOT temporal variability, the AERONET network is the better candidate, due to its 

166 frequent temporal sampling during daylight hours, unlike the snapshots from MISR and MODIS. 

167 However, when the spatial and temporal variability is examined simultaneously, the MISR and MODIS 

168 AOT retrievals can also provide useful information about space-time variability. For simplicity, the 

169 spatial and temporal analysis is presented here using the MISR and MODIS data, but the conclusions 

170 about the temporal variability are consistent with those obtained using the AERONET observations 

171 (results not shown). The temporal component of the analysis is useful for identifying the timescales over 
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172 which the AOT data can be integrated into relatively contiguous maps without introducing errors due to 

173 correlations in the temporal variability of the AOT. 

174 The spatiotemporal correlation analysis is performed using variogram analysis (e.g., Chiles and Delfiner, 

175 1999). For all pairs of AOT data from a given instrument (e.g. MISR), the raw variogram is evaluated as: 


1 


176 where z are the AOT observations at locations jc, and Xj and times U and tj, h x is the spatial separation 

177 distance between the two observation locations, and h t is the temporal lag in days between the 

178 observations. h x is calculated as the great circle distance between the locations jc, and xj 


( 2 ) 


K ( x i , x j) - r cos 1 sin (j>. sin (f>, + cos <j> t cos <f>, cos(d?. - 0 . ) 


( 3 ) 


179 where ( <f>u d t ) are the longitude and latitude of location x t and r LV □ WKH □ (DUWKfV □ SHdljlfiLX V n □ 

180 presented here, a raw variogram is created for each repeat cycle of MODIS and MISR (i.e. each available 

181 16 day period in 2001). 


1 82 Once the raw variograms are obtained, the variability can be visualized by binning the variances y into 

1 83 preset ranges of separation distances (h x ) and time lags (h t ). The binned version of the raw variogram is 

1 84 referred to as the experimental variogram. If the temporal correlation of the observations across multiple 

1 85 days is negligible, the experimental variogram can be presented as a function of spatial lag only, and a 

1 86 theoretical model can be selected to represent the observed spatial-only variability. In the analyses 

1 87 presented here, an exponential model was found to represent the spatial correlation of the AOT data well: 
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2 2 - ( 4 ) 

188 where a 1 (=cr 2 Co \ ) represents the variance of observed AOT at large separation distances (i.e., for 

1 89 uncorrelated observations) and l is the range parameter. The correlation length beyond which correlation 

190 between points becomes negligible is defined as approximately 31 [e.g., Chiles and Delfiner, 1999]. cr 2 is 

191 the nugget, representing both the measurement error and the small-scale variability at distances smaller 

192 than those resolved by available observations, whereas cr 6 2 represents the variance of the portion of the 

193 AOT variability that is spatially correlated. These parameters are optimized using a least squares fit to the 

194 spatial raw variogram. Conceptually, a higher variance is indicative of greater overall variability, and a 

195 shorter correlation length indicates greater spatial variability at smaller scales. 

196 3.3 Geostatistical Data Fusion Approach 

197 Universal kriging (e.g., Chiles and Delfiner, 1999), a geostatistical data fusion approach, makes it 

198 possible to fuse auxiliary variables with full spatial coverage (e.g. MISR and MODIS AOT) to improve 

199 the interpolation of a primary dataset with observations at a finite number of locations (e.g. AERONET 

200 AOT). The auxiliary variables fill a role analogous to regressors in multiple linear regression, but within a 

201 framework that accounts for the spatial correlation of the estimated field, and can reproduce observed 

202 AERONET AOT measurements exactly at sampling locations. 

203 The objective is to estimate the AOT distribution (s) at m locations and times (typically defined on a 

204 regular grid), given the AERONET AOT measurements at n locations and times, where s (m x 1) is 
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modeled as the sum of a deterministic but unknown component X, 



220 AOT observations. At the estimation locations, this component represents the predicted residuals between 

221 the true AOT and the weighted MISR and MODIS AOT at those locations/times. The covariance of these 

222 residuals is described using a matrix Q, where the covariance function is defined based on the variogram 

223 analysis (Equation 4), such that the covariance between two points x, and Xj is defined as: 

2 

2 2 

( 7 ) 

2 

In this case, the variogram analysis is carried out on the detrended AERONET AOT data (z ± X s 
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235 where Q zz is the n x n spatial covariance matrix defined between AERONET observation locations 

236 based on Equation 7, X z (n x p) is the trend term defined at the measurement locations based on Equation 

237 6, Q zs (n x m) represents the covariance evaluated between the measurement and the estimation locations 

238 again based on Equation 7, and X s (m x p) is the model of the trend defined at the estimation locations 

239 again based on Equation 6.When multiple time periods are used in the analysis, as is true for the test cases 

240 described in Section 3.4, the correlation between time periods is assumed to be zero. The system of 
equations is solved for M, a p x m matrix of Lagrange multipliers, and 
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248 where the diagonal elements of 0 represent the uncertainty of the individual drift parameters, and the 

249 off-diagonal terms represent the estimated covariance of the errors associated with these estimates. Recall 

250 that the drift coefficients are the weights assigned to the MISR and MODIS AOT datasets. These weights 

25 1 remain constant over the domain of analysis. As a consequence, the relationship between the true (as 

252 represented by AERONET) AOT, and MISR and MODIS AOT, is implicitly assumed to remain constant 

253 within an examined region. This is one of the reasons for which the test cases examined in Section 3.4 are 

254 conducted regionally, because the relationship between MISR, MODIS, and AERONET cannot 

255 necessarily be assumed to remain constant throughout the continental United States. 

256 Ordinary Kriging, a simple geostatistical interpolation technique, is used for comparison to the 

257 Universal Kriging estimates in the presented analyses. Ordinary kriging is one of the most commonly 

258 used techniques in geostatistical gap filling, but it lacks the advantage of using information from multiple 

259 sensors. In the ordinary kriging approach, X s = [l...l] r , and the covariance is derived using a variogram of 

260 the AERONET observations without detrending. The other equations remain unchanged. Past 

261 applications of ordinary kriging in aerosol science have been limited to estimation of aerosol species over 

262 different regions [Zapletal, 2001 ; Delalieux et al., 2006], and have not been aimed at comparison with 

263 other estimation techniques. In this work, because the true AOT distribution is unknown, the ordinary 

264 kriging estimates are used as a baseline for evaluating the estimates from universal kriging. By comparing 

265 the two kriging estimates, we identify the effect of using additional satellite observations on both the 

266 AOT estimates and the uncertainty associated with those estimates. 

267 3.4 Test Cases 

268 The correlation analysis (Section 3.1) is carried out for three regions over the contiguous United States for 

269 selected periods in 2001 . Recognizing that aerosol distributions can be both site and season specific, the 
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270 United States are divided into three regions (Western, Central, and Eastern), as illustrated in Figure 1. In 

271 the Western region, we expect dust to be dominant, along with biomass burning during the summer and 

272 autumn months. Biogenic aerosols often dominate the southeast, especially in summer, where biomass 

273 burning may also be important in some seasons. Four seasons are considered: Winter (DJF), Spring 

274 (MAM), Summer (JJA) and Fall (SON). Following previous studies ( Kahn et al. [2005a]; Liu et al. 

275 [2004a]; Abdou et al. [2005]), the correlation analysis was carried out at a daily scale. 

276 The spatio-temporal analysis is carried out using the MISR and MODIS datasets. The analysis is 

277 performed at the native resolution of MISR (i.e., 17.6 km) and MODIS (i.e., 10 km) AOT for each 16-day 

278 repeat cycle of the Terra satellite in 2001 . By doing this analysis for each 16-day repeat cycle, the 

279 seasonal changes in spatio-temporal variability of AOT can be assessed, as well as how these changes 

280 relates to the periodic changes in the underlying AOT processes over the continental US. 

281 Finally, the geostatistical data fusion analysis is presented using two test cases, the first being over the 

282 Eastern US during autumn, and the second over the Western US during summer (Figure 2). Table 1 

283 outlines the details of these two test cases. For these test cases, the study region is broken up into 0.2° x 

284 0.2° grid cells at which the AOT estimates are obtained. The MISR and MODIS observations used for 

285 this analysis are averages of all the MISR and MODIS AOT observations falling within a given 0.2° x 

286 0.2° grid cell. This particular estimation resolution is chosen to show the flexibility of the universal 

287 kriging approach in estimating AOT at very fine resolutions, but in general could be performed at coarser 

288 estimation scales as well. 

289 As will be shown in Section 4.2, there is little significant temporal correlation in the day-to-day variability 

290 in the MISR and MODIS AOT within a 7-day period. As a result, the geostatistical data fusion is 

291 performed in one-week increments, using weekly-averaged AOT data from AERONET, MISR and 
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292 MODIS. These averaged AOT data are used to obtain estimates of the average spatial distribution of AOT 

293 over those 7-day periods. For each 7-day period, AERONET sites that have AOT data for at least three of 

294 the seven days, and which have overlapping MISR and MODIS data, are used in the analysis. As a result, 

295 for the Eastern Test Case during autumn, two to ten sites are used during the various weeks, whereas for 

296 the Western Test Case during summer, two to eight sites are used in each week. Figure 2 shows the 

297 locations of all the sites used in the test cases. It should be noted here that there may be cases in which 

298 significant temporal correlation may exist (e.g. near sources), such as where shorter time-scale variations 

299 are predominant, or where strong gradients occur in transported aerosol from far sources. In such cases 

300 the data fusion approach should be applied with caution. 

301 The AOT estimates obtained from universal kriging (henceforth denoted as AOTuk) are compared with 

302 AOT estimates obtained from ordinary kriging (henceforth denoted as AOT 0 k). Cross-validation is used 

303 to compare the two estimates. In this approach, individual AERONET 7-day observations at a given site 

304 are sequentially eliminated from the analysis, and estimates at these locations and times are obtained 

305 using the remaining AERONET observations, and, for AOTuk, using available MISR and MODIS data as 

306 well. Because AERONET measurements have traditionally been used for validating satellite observations 

307 of MISR and MODIS [Kahn et al., 2005a; Renter et al., 2005; Yu et al., 2006], the withheld AERONET 

308 observations are used to evaluate the relative precision and accuracy of the AOT 0 k and AOTuk estimates. 

309 The evaluation of AOT 0 k and AOTuk estimates is carried out using three metrics. First, the root mean 

310 square error (RMSE) is calculated between the estimated AOT and the AERONET observations. Second, 

311 the magnitude of the predicted kriging uncertainties is evaluated by calculating the root mean square 

312 prediction error (RMSPE) of the kriging uncertainties (Equation 10). Third, the accuracy of these 

313 predicted uncertainties is evaluated by verifying the percent of true AERONET AOT observations that 
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314 fall within two standard deviations of the estimated AOT, where the standard deviations are those 

315 predicted by the kriging analyses (Equation 10). This third metric is less sensitive to extreme outliers, 

316 and, in an ideal scenario, 95% of the true AOT should fall within this interval. Values significantly below 

317 95% would indicate an underestimation of the true uncertainties, while values substantially above 95% 

318 indicate overly conservative estimates. All three metrics are calculated across the entire season for both 

319 test cases. 

320 Overall, the two examined test cases are designed to (i) demonstrate the versatility of the universal 

321 kriging technique in estimating AOT over different regions and across seasons, (ii) evaluate the 

322 improvement of universal kriging estimates over ordinary kriging estimates (or simply the AOT fields 

323 observed by MISR or MODIS individually) as a function of the strength of the relationship between 

324 MISR, MODIS, and AERONET AOT. 

325 4 RESULTS AND DISCUSSION 

326 4.1 Comparison of MISR, MODIS and AERONET datasets 

327 The results of the correlation analysis are presented in Table 2, and reveal that MISR data have a stronger 

328 correlation to AERONET data (p=0.47 to 0.92) than do MODIS data (p=0.09 to 0.61) across seasons in 

329 the Western and Central regions. Given that the correlation coefficient is an indication of the degree of 

330 linear covariability between the datasets, this implies that MISR is better able to explain the variability in 

331 the AERONET AOT than MODIS over these regions. Note that for MODIS, the "standard" product was 

332 used in this paper, and it is plausible that using the more recent "Deep Blue" product, which was not 

333 available at the time of analysis, would have shown better agreement with the AERONET AOT, at least 

334 over bright surfaces. On the other hand, the particle properties used in the MODIS standard AOT retrieval 
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335 over land are assumed based on AERONET values [Levy et al., 2007], whereas for MISR, particle 

336 properties are retrieved along with AOT as part of a self-consistent process [Martonchik et al., 2009]. 

337 The weak correlation in the Western region for both instruments is primarily due to low AOT values in 

338 this region, which are near the lower limit of retrieval sensitivity for MISR and MODIS. Liu et al. [2004a; 

339 2004b] points out that low values of AOT, as well as coarse-particle dominated scenarios, may produce 

340 poor correlations with AERONET AOT. This does not necessarily indicate poor MISR or MODIS 

341 performance; rather, at very low AOT values, the correlation coefficients are not informative because the 

342 uncertainty associated with the satellite retrievals is large compared to the magnitude of the AOT itself. 

343 The high surface albedo in the Western sector and the frequent atmospheric loading with non-spherical 

344 mineral dust are additional obstacles to obtaining good satellite retrievals of AOT over this region. 

345 In the Eastern region, the MISR (p=0.52 to 0.86) and MODIS (p=0.70 to 0.87) data show comparable 

346 correlations to AERONET across seasons, and are able to capture the AERONET AOT variability better 

347 than across the other two examined regions. The year-round correlation coefficients (p = 0.78 for MISR 

348 and p=0.84 for MODIS) are similar to values that have been reported previously for continental sites in 

349 this region [Chu et al., 2002; Liu et al., 2004a; Kahn et al., 2005a]. 

350 Overall, results from this analysis are consistent with previous findings that indicate that differences 

351 between MISR and MODIS AOT relative to AERONET AOT are caused by site-specific effects and 

352 aerosol-size-distribution effects. 

353 This initial analysis indicates that the information provided by MISR and MODIS with regard to the AOT 

354 distribution as measured by the AERONET network varies regionally and seasonally throughout the 

355 continental United States. Based on these results, it is expected that universal kriging analysis should 
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356 outperform ordinary kriging in the Eastern region, where the correlations between the AERONET data 

357 and the MISR and MODIS data are strong. In other regions, the additional information provided by 

358 MODIS and MISR is less significant, and the universal kriging and ordinary kriging analyses are 

359 expected to be more similarly to one another. 

360 4.2 Spatio-Temporal Variability Analysis 

361 The spatio-temporal variability analysis is performed for MISR and MODIS AOT for each 16-day repeat 

362 cycle in 2001 for the Terra satellite. For each repeat cycle, a spatio-temporal experimental variogram is 

363 obtained. Example variograms are presented for MISR and MODIS in Figures 3a and 3b, for April 1 1 to 

364 26, 2001. These variograms represent the expected variance of pairs of MISR (Figure 3a) or MODIS 

365 (Figure 3b) observations, separated by a given distance in space and time lag. 

366 Figures 3a and 3b do not exhibit any noticeable temporal correlation in the day-to-day variability of the 

367 AOT distribution for time lags up to 7 days. The temporal lag (shown on the vertical axis) in Figures 3a 

368 and 3b represents the number of days between the times when two observations are recorded. A temporal 

369 lag of one day could therefore represent, for example, the expected variance between observations taken 

370 on days 14 and 15, or on days 1 and 2 or the repeat cycle. The lack of temporal correlation indicates that 

371 the coherent temporal variability in the AOT takes place either at sub-diumal scales that cannot be 

372 captured by the examined remote sensing data products, and/or at longer time scales, potentially 

373 representative of seasonal variability. 

374 Hence, the results of this analysis indicate that temporal correlation is not significant at for time lags up 

375 to 7 days, and therefore that MISR and MODIS data taken over a week can be integrated into a single 

376 map. In other words, multiple days of data can be used concurrently to inform the data fusion analysis. 
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377 Note that the lack of temporal correlation does not necessarily imply a lack of temporal variability, but 

378 simply that the observed AOT is not correlated from day to day. 

379 On the other hand, Figures 3a and 3b reveal that the MISR and MODIS AOT data do exhibit strong 

380 spatial correlation, as evidenced by the fact that the variance increases as the spatial separation distance 

381 increases. This is more clearly visible in Figures 3c and d, where all data from the April 1 1 to 26 repeat 

382 cycle are examined in a single spatial variogram. These figures display both the experimental and the 

383 fitted theoretical spatial variograms for MISR and MODIS. These two figures indicate that the correlation 

384 length of AOT data (i.e., the lag distance at which the semivariance reaches an asymptote) appears to be 

385 approximately 900km for both instruments for the examined time period, indicating that observations 

386 separated by longer distances are essentially independent. 

387 Figures 3e and 3f present the parameters of the fitted theoretical spatial variograms for each of the 25 

388 Terra repeat cycles in 2001 . This analysis shows that the spatial correlation of the MISR and MODIS 

389 AOT data are quite consistent with one another (blue lines in Figures 3e and f). The correlation lengths 

390 vary significantly throughout the year, ranging from 500km to 1 500km, with higher values prevalent 

391 during the summer months. On the other hand, the total amount of variability (i.e. variance) of the 

392 MODIS AOT is always significantly higher than that of MISR. During the winter months, both MISR and 

393 MODIS show shorter correlation lengths and increased variance, representative of a more heterogeneous 

394 distribution of aerosols. In general, the long-range transport of dust in late spring and summer, and smoke 

395 from summer through early autumn, are likely to contribute to the longer correlation lengths during the 

396 summer months, whereas local aerosol sources explain the smaller-scale variability observed during other 

397 seasons. 
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398 Seasonal changes in the spatial variability of AOT will impact the uncertainty estimates obtained from 

399 universal kriging. During the summer months, due to the longer correlation lengths and smaller variance, 

400 the AOT estimates will have lower uncertainty, while, conversely, during the winter months, we can 

401 expect higher estimation uncertainties. 

402 One interesting conclusion from Figure 3 is that the MODIS AOT variance is higher than that for MISR, 

403 across all seasons. This is due in part to the fact that the more frequent and finer scale MODIS sampling 

404 captures more small-scale AOT variability than MISR. The second reason for this higher variance is that, 

405 due to its exclusively near-nadir viewing geometry, MODIS has a greater sensitivity to variability in 

406 surface brightness on small spatial scales, which in turn introduces some additional variability into the 

407 MODIS AOT retrievals. Neither of these features hinders the application of the universal kriging 

408 approach presented in this work. However, it has implications for researchers pursuing assimilation of 

409 MISR and MODIS radiance data, or looking to improve the retrieval algorithms of these two sensors. 

410 4.3 Data Fusion Results 

411 Figures 4 and 5 show the estimated AOT field for one week of each of the case studies described in Table 

412 1 . The Eastern test case demonstrates that the universal kriging AOT estimates are better than the 

413 ordinary kriging estimates when MISR and MODIS are significantly correlated with the AERONET AOT 

414 observations. The associated uncertainties for the AOTmc estimates are significantly lower. Cross- 

415 validation at the AERONET locations confirms that the AOTmc estimates are more realistic than the 

416 AOT 0 k estimates, as shown for one of the 7-day periods in Figure 6 (see Figures SI and S2 in the 

417 Supplementary material for the entire season). Overall, for this test case, the RMSE for AOTmc is 0.053, 

418 which is lower than that of AOT 0 k (0.067) and each of the individual satellite datasets (0.054 for MISR 

419 and 0.056 for MODIS). The true AOT falls within the 2 standard deviations of both the kriging estimates 
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420 for 93% of AERONET observations, but the RMSPE of AOTuk (RMSPE = 0.035) is significantly lower 

421 than that of AOT 0 k (RMSPE = 0.069). These results confirm that, when strong correlation exists between 

422 multiple datasets, the universal kriging approach can be used to obtain better predictions with smaller 

423 uncertainties relative to estimates based on measurements from a single sensor. This is evident not only 

424 from the reduction in uncertainty, but also from the lower RMSE and RMSPE values of AOTuk relative 

425 to AOTqk. 

426 The Western test case demonstrates that the universal kriging estimates are comparable to the ordinary 

427 kriging estimates in regions where the correlation with MISR and MODIS is low. The predicted 

428 uncertainty (Figure 5b and d) is similar for the two methods. This is consistent with our findings from the 

429 correlation analysis,, because MISR and MODIS are not strongly correlated with the AERONET AOT in 

430 this region (Table 2), and are therefore unable to capture the AERONET AOT variability. Cross- 

431 validation results shown in Figure 7 confirm that the two approaches provide similar estimates with high 

432 uncertainty (see Figure S3 and S4 in the Supplementary material for the entire season). The AERONET, 

433 MISR and MODIS AOT values are also plotted in Figure 7, and demonstrate that, for the examined case, 

434 both ordinary and universal kriging seem to do better than just using individual MISR and MODIS 

435 datasets. This is further validated by the metrics calculated for the entire season. The RMSE for both 

436 AOTuk and AOT 0 k is 0.047 and 0.048, respectively, which is lower than the RMSE of 0.082 for MISR 

437 and 0.26 for MODIS. The true AERONET AOT fall within 2 standard deviations of AOTuk and AOT 0 k 

438 estimates for 98% of available observations. Finally, the RMSPEs are similar for the ordinary (RMSPE = 

439 0.066) and universal (RMSPE = 0.061) kriging approaches, reaffirming their similarity to one another for 

440 this test case. 
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441 In addition to predicting AOT, the universal kriging approach can be used to quantity which of the 

442 satellite observations has more influence on the estimation procedure, by looking at the drift coefficient 

443 ( O) values and their uncertainties ( <Tq), as shown in Table 3. A coefficient of variation ( a below 

444 0.5 implies a statistically significant contribution to the trend at the 2 level. For the Eastern Test Case, 

445 the MODIS AOT observations have a more significant drift coefficient (<Tq/ 0.18) than the MISR 

446 data, and these latter data are therefore used primarily to adjust the spatial pattern in MODIS AOT to 

447 more closely resemble the AERONET AOT observations. Conversely, for the Western Test Case, the 

448 MISR AOT observations seem to be a significant predictor of AERONET AOT measurements ( a q/ ® = 

449 0.43). In addition, the drift coefficient values for the constant term ( O ) are not significantly different 

450 from zero for either examined case, indicating an absence of any systematic offset between the AOT 

451 predicted by the weighted combination of MISR and MODIS, and the AOT observed by AERONET. 

452 Finally, although this analysis used both MISR and MODIS in the data fusion process, one could easily 

453 use either MISR or MODIS individually, or some other combined AOT product(s). The approach 

454 presented combines the best available information from all available sensors to identify the optimal 

455 weighted combination to represent the AOT distribution. 

456 5 CONCLUSIONS 

457 A geostatistical data fusion technique is implemented for combining remote-sensing and ground-based 

458 observations of AOT. Results show that adopting the universal kriging approach based on the 

459 combination of MISR, MODIS and AERONET enables better estimation of AOT with reduced 

460 uncertainties, relative to estimates based on observations from a single instrument. 
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461 All three examined datasets were found to display strong spatial correlation in their measured AOT 

462 distributions. Although the total degree of AOT variability differed between MISR and MODIS, the 

463 spatial scales of this variability were similar for these instruments. The day-to-day temporal correlation in 

464 MISR or MODIS AOT observations was found to be minimal, due at least in part to limited temporal 

465 sampling, making it possible to integrate such observations over multiple days to better infer the spatial 

466 distribution of AOT. 

467 As an increasing number of remote sensing observations become available, data fusion approaches such 

468 as the one presented here may hold the key to furthering our understanding of atmospheric aerosols. 

469 Although differences between instruments are always present, the approach implemented here takes 

470 advantage of their complementary features by combining the datasets in a manner that is statistically 

471 robust. The approach relies on the availability of auxiliary variables (MISR and MODIS AOT, in the 

472 presented analysis) at all locations where the AOT is to be estimated, and assumed that the relationship 

473 between these auxiliary variables and the primary observations (AERONET AOT, in the presented 

474 analysis) remains constant throughout the examined region. 

475 Finally, this study reinforces the complementary value of remote-sensing and ground-based observations 

476 of AOT. Long-term monitoring of aerosol distributions is possible via remote sensing measurements, and 

477 these can be used to capture the spatio-temporal distribution of aerosols. Expected refinements in retrieval 

478 algorithms and sensor capabilities will improve the accuracy of the retrieved AOT further. By fusing 

479 these measurements with ground-based observations using techniques such as the one presented here, it 

480 will be possible to obtain reliable long-term estimates of AOT at national and even global scales. 
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638 TABLES 

639 

640 Table 7- Test case specifications 


Test 

Case 

Time 

Period 

Spatial 

Extent 

Estimation Resolution 

Number of 
AERONET 
locations 

Spatial 

Temporal 

Eastern 

Fall 

70°-85°W 

25°-50°N 

0.2° x 0.2° 

Average 
over a 7 day 
period 

10 

Western 

Summer 

105°-120° W 
25°-50° N 

0.2° x 0.2° 

Average 
over a 7 day 
period 

8 


641 

642 Table 2 ± Correlation coefficients between AERONET measurements and MISR and MODIS 

643 observations classified by region (Figure 1), and season for the year 2001 . Low correlation coefficient (0- 

644 0.5) cases are shaded in dark gray, medium correlation coefficient (0.5-0.75) cases are shaded in light 

645 gray, and high correlation coefficient cases (0.75-1 .00) are in bold. The lowest correlations occur in the 

646 west, where bright surfaces and mixtures of spherical particles and non-spherical dust dominate, and in 

647 the winter months, when total-column AOT tends to be low, and AOT is near the sensitivity limit of the 

648 satellite instruments. 



Winter (DJF) 

Spring (MAM) 

Summer (JJA) 

Fall (SON) 

All Months 

MISR 

MODIS 

MISR 

MODIS 

MISR 

MODIS 

MISR 

MODIS 

MISR 

MODIS 

Western 

0.49 

0.09 

0.67 

0.29 

0.47 

0.32 

0.59 

0.09 

0.63 

0.30 

Central 


0.51 

0.68 

0.61 

0.73 

0.58 

0.67 

0.43 


0.59 

Eastern 


0.70 

0.52 









649 

650 Table 3 - Drift coefficient values, their associated uncertainties and the coefficient of variation for the two 

651 test cases from the universal kriging model. 


Test 

Case 

Constant 

MISR 

MODIS 

O 

0 


o-q/^o 

O 

i 



O 

2 



Eastern 

0.010 

0.014 

1.4 

0.017 

0.16 

9.42 

0.68 

0.12 

0.18 

Western 

0.027 

0.024 

0.88 

0.37 

0.16 

0.43 

0.040 

0.070 

1.75 
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653 FIGURES 

654 

655 Figure Captions 

656 

657 Figure 1 - Location of AERONET sites used in the correlation analysis. The examined regions are 

658 outlined in blue. Closed circles represent observation locations also used in the universal kriging test 

659 cases. 

660 

661 Figure 2 ± AERONET sites used for Western and Eastern test case. The AERONET sites in the Western 

662 Region are Missoula (MIS), Rimrock (RIM), BSRN-BAO-Boulder (BSR), Railroad Valley (RAI), Rogers 

663 Dry Lake (ROG), La Jolla (LAJ), Maricopa (MAR) and Sevilleta (SEV); the AERONET sites in the 

664 Eastern Region are Rochester (ROC), Cartel (CAR), Harvard Forest (HVF), GISS (GIS), Philadelphia 

665 (PHI), MD Science Centre (MDS), GSFC (GSF), Big Meadows (BIG), Wallops (WAP) and Cove (COV). 

666 

667 Figure 3 - Variograms of the spatial and temporal variability of AOT from MISR and MODIS. Panels (a) 

668 and (b) represent the spatial and temporal variograms of AOT over a 16 day period from April 1 1 to April 

669 26, 2001 for MISR and MODIS respectively. The color bar indicates the semi-variance. Panels (c) and (d) 

670 represent the spatial variogram over the same period from MISR and MODIS, respectively, using all data 

671 over the 16-day period. Panels (e) and (f) show the correlation length ( 31) and variance of AOT for all 16- 

672 day periods in 2001 for MISR and MODIS, respectively. Note that MODIS had no data for one repeat 

673 cycle in June. 

674 

675 Figure 4 - Comparison of AOT 0 k with AOTmc for Eastern Test Case for one period from October 29 to 

676 November 4. The black asterisks indicate the locations of the AERONET sites. The white gaps indicate 

677 the 7-day satellite coverage mask that is imposed on both universal and ordinary kriging for ease of 

678 comparison, (a) Best estimates obtained from ordinary kriging. (b) Uncertainty associated with the 

679 ordinary kriging estimates, (c) Best estimates obtained from universal kriging. (d) Uncertainty associated 

680 with the universal kriging estimates. 

681 

682 Figure 5 - Comparison of AOT 0 k with AOTmc for Western Test Case for a 7-day period from July 21 to 

683 July 27. The black asterisks indicate the locations of the AERONET sites. The white gaps indicate the 7- 

684 day satellite coverage mask that is imposed on both universal and ordinary kriging for ease of 

685 comparison, (a) Best estimates obtained from ordinary kriging. (b) Uncertainty associated with the 

686 ordinary kriging estimates, (c) Best estimates obtained from universal kriging. (d) Uncertainty associated 

687 with the universal kriging estimates. 

688 

689 Figure 6 ± Cross validation results for October 29 ± November 4 for Eastern test case. Error bars 
UHSUHVHQWn □ 
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694 SUPPLEMENTARY MATERIAL 

695 

696 Figure Captions 

697 

698 Figure SI ± Cross-validation results for September 3 to October 14, 2001 for Eastern test case. Error bars 
UHSUHVHQW n □ 
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